unit MatrixOperations;

{
  Матричные операции
  (с) Keeper, 2011
}

{
  PasCalc
  Version 1.2 2011-05-26

  Copyright (C) 2011 KEEPER <yar.keeper@gmail.com>

  This file is part of the PasCalc Calculator

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.

  See the file LICENSE.TXT, included in this distribution,
  for details about the copyright.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
}

interface

uses
  Math, SysUtils;

type

  TVector = array of extended; //строка матрицы

  TMatrix = array of TVector;  //сама матрица

{ Объявление исключений }

  { Требуется ненулевая матрица }
  ENotNullMatrixExcepted = class(Exception)
    constructor Create;
  end;

  { Матрица имеет неверный формат для применения метода Гаусса }
  EIllegalMatrixFormat = class(Exception)
    constructor Create;
  end;

  { Система не имеет решения }
  ENotSolution = class(Exception)
    constructor Create;
  end;

{ Объявление процедур и функций }

{
  Решаем СЛАУ методом Гаусса. SourceMatrix - матрица вида n x n+1
  n+1 столбец - столбец свободных членов СЛАУ
  Функция возвращает вектор - столбец решений
}
function GausseMethod(const SourceMatrix: TMatrix): TVector;

{ Преобразует матрицу SourceMatrix в треугольную }
procedure ConvertMatrixToTriangle(var SourceMatrix: TMatrix);

{ Возвращает копию матрицы SourceMatrix }
function CopyMatrix(const SourceMatrix: TMatrix): TMatrix;

implementation

function CopyMatrix(const SourceMatrix: TMatrix): TMatrix;
var
  i, j: integer;
begin
  SetLength(Result, Length(SourceMatrix), Length(SourceMatrix[0]));
  for i := 0 to high(SourceMatrix) do
    for j := 0 to high(SourceMatrix[i]) do
      Result[i, j] := SourceMatrix[i, j];
end;

{ Получаем вектор с решением СЛАУ }
function GetSolution(const TriangleMatrix: TMatrix): TVector;
var
  i, j: integer;
  RowCount, ColCount: integer;
  sum: extended;
begin
  RowCount := Length(TriangleMatrix);
  ColCount := Length(TriangleMatrix[0]);
  if (TriangleMatrix[RowCount - 1, ColCount - 1] = 0) and
    (TriangleMatrix[RowCount - 1, ColCount - 2] = 0) then
      raise ENotSolution.Create;
  SetLength(Result, RowCount);
  result[RowCount - 1] := TriangleMatrix[RowCount - 1, ColCount - 1]/
    TriangleMatrix[RowCount - 1, ColCount - 2];
  for i := RowCount - 2 downto 0 do
  begin
    sum := TriangleMatrix[i, ColCount - 1];
    j := ColCount - 2;
    while j > i do
    begin
      sum := sum - (Result[j] * TriangleMatrix[i, j]);
      dec(j);
    end;
    Result[i] := sum/TriangleMatrix[i,j];
  end;
end;

function GausseMethod(const SourceMatrix: TMatrix): TVector;
var
  WorkMatr: TMatrix; //рабочая матрица над которой производятся преобразования
begin
  if SourceMatrix = nil then
    raise ENotNullMatrixExcepted.Create;
  //Проверяем возможность решения СЛАУ
  if Length(SourceMatrix) <> Length(SourceMatrix[0]) - 1 then
    raise EIllegalMatrixFormat.Create;
  //Создаём рабочую матрицу чтобы не портить исходные данные
  WorkMatr := CopyMatrix(SourceMatrix);
  //Приводим рабочую матрицу к треугольному виду
  ConvertMatrixToTriangle(WorkMatr);
  //Получение решения
  result := GetSolution(WorkMatr);
end;

{ Преобразует матрицу в треугольную }
procedure ConvertMatrixToTriangle(var SourceMatrix: TMatrix);
var
  i, j, k: integer;
  coef: extended;
begin
  for i := 0 to High(SourceMatrix) - 1 do
    for j := i + 1 to high(SourceMatrix) do
    begin
      coef := -(SourceMatrix[j, i]/SourceMatrix[i, i]);
      for k := i to high(SourceMatrix[j]) do
        SourceMatrix[j, k] := SourceMatrix[j, k] + (SourceMatrix[i, k] * coef);
    end;
end;

{ ENotNullMatrixExcepted }

constructor ENotNullMatrixExcepted.Create;
begin
  inherited Create('Matrix must be not null!');
end;

{ EIllegalMatrixFormat }

constructor EIllegalMatrixFormat.Create;
begin
  inherited Create('Matrix has illegal format. Matrix must have n rows and n+1 cols');
end;

{ ENotSolution }

constructor ENotSolution.Create;
begin
  inherited Create('System of equations has no solutions');
end;

end.
